


###——————————————————————背景知识——————————————————————————
#生存时间通常不呈正态分布；普通的线性回归不能有效利用删失数据；
#目前生存时间最常用的模型是Cox比例风险回归模型。它不要求估计资料的生存分布类型。
#在R中，strata()用于识别方程右边的分层分类变量。如strata(age);class

###基本概念
#生存函数：S(t,X)，具有变量X的观察对象，经历t个时间段后，仍存活的可能性。死亡函数F(t,X)=1-S(t,X)。
#死亡密度函数：f(t,X)，具有变量X的观察对象，在时刻t的瞬时死亡率。f(t,X)=F'(t,X)。
#风险函数：h(t,X)，代表在t时刻仍存活的个体（具有变量X），在t时刻的瞬时死亡率。h(t,X)=f(t,X)/S(t,X).

###模型概况
#Cox回归模型的基本形式——比例风险率模型( proportional hazards model,PH Model)：ln[h(t,X)/h0(t)]=lnRR=β1X1+β2X2+β3X3+β4X4+……βiXi+……
###Cox回归应用的条件
#比例风险假定：各危险因素的作用不随时间的变化而变化，即h(t,X)/h0(t) 不随时间的变化而变化。
#对数线性假定：模型中的自变量应与风险比的对数呈线性关系。

##PH假定判断：满足PH假定的资料才能进行Cox回归分析。
#1.cox-km图法比较：通过观察某变量分组的KM生存曲线，如果交叉，提示不满足PH假定。交叉时应该选用landmark分析。
#2.lsl或lml图：绘制某变量分组的ln[-ln(生存率)]对生存时间t的图，曲线大致平行或等距，则各协变量均满足或近似满足PH假定，可直接应用基本Cox模型。此过程，如果有连续变量（如年龄），则需要转变为二分量变量。

##R语言中的Cox回归模型诊断
#比例风险假定：Schoenfeld残差
#模型影响点（异常值）识别：Deviance残差
#比例风险的对数值与协变量之间的非线性关系识别；Martingale残差
#如果违反比例风险假设可以通过以下方式解决：模型中添加协变量与时间的交互相应；分层分析。

##偏回归系数β的解读
#偏回归系数βi的流行病学含义是:在其他协变量不变的情况下，协变量Xi每增加一个测定单位时所引起的相对危险度的自然对数的改变量。
#若βi>0,RRi>1，则Xi取值越大时，h(t,X)的值越大，此时Xi为危险因素。若βi=0,RRi=1，则Xi取值对h(t,X)没有影响，此时Xi为无关因素。若βi<0,RRi<1，则Xi取值越大时，h(t,X)的值越小，此时Xi为保护因素。即β>0危险，β<0保护。
#exp(βi)即代表协变量Xi每增加一个测定单位时所引起的相对危险度的改变量。
#通常，对于某个观察对象，Cox回归方程还可获得：PI。右边的β1X1+β2X2+β3X3+β4X4+……βiXi+……称为预后指数（prognostic index,PI）。PI越大，预后越差。预后指数可直观比较不同患者的预后，或者用于低中高危分组等，对于指导治疗、提高生存率有重要作用。估计生存率。S(t)=[S0(t)]的exp(PI)次方，S0(t)用Breslow估计法。比如，肿瘤分级I级，肿瘤<3.0cm，未复发，59个月的生存率为25.6%。

##变量筛选：
#类似于多重线性回归与logistic回归，有Forward法、Backward法和逐步回归法（Stepwise法）。通常的，α引入=0.05,α剔除=0.10?

##参数估计：
#Cox回归模型中的偏回归系数的估计通常借助部分似然理论（partial likelihood）,因为第2部分只与β有关，这样则不需要确定基准函数h0(t)的值；另外部分似然只与生存时间的排序而不是大小有关，这样生存时间的单调变换不会影响β的估计值。通常利用Newton-Raphson迭代法求得β。
#通常需要获得β、SE(β)、Wald-X2、P、RR、RR95%CI、标准化β等值。其中，标准化β可以比较协变量的影响大小。

##模型校正：
#在建模过程中，需要仔细考虑不同协变量的共线性问题。比如，肿瘤患者预后与分级、年龄、肿瘤大小、是否复发等协变量相关。按医学常识，年龄在肿瘤分级中很可能是均衡分布的，但肿瘤大小、是否复发等在肿瘤分级中的分布显然是不均衡的，因为肿瘤越大常常预示着更差的肿瘤分级，而且复发概率也与肿瘤分级有很大的相关性，即协变量音的多重共线性问题。此时，如果不进行校正，最后由于肿瘤大小和是否复发的“掩盖”下，肿瘤分级会被认为对肿瘤患者预后没有影响，这显然违背了医学的常识。因此，在Cox回归建模过程中，需要仔细考虑不同协变量的共线性问题。在R中，可通过(X1 + X2 + X3) * X4的形式来处理可能存在多重共线性的变量以进行Cox回归模型的校正。

##假设检验：
#回归模型的假设检验通常采用似然比检验、Wald检验和Score检验，其检验统计量均服从X2分布，其自由度为模型中待检验的自变量个数。一般说来，Cox回归系数的估计和模型的假设检验计算量较大，通常需利用计算机来完成相应的计算。


###————————————————————工作环境————————————————————————————
setwd("E:/RCloud/R learning metirial/sample document/COX regression")
library(survival);library(survminer);library(car)



###——————基于lung这一survival自带data进行示例————————————

data(lung)#NCCTG Lung Cancer Data
set.seed(2018);lung[sample(1:228,20),]#

res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data = lung)
summary(res.cox)

#进行PH假定：如结果，三个变量的P值都大于0.05，说明每个变量均满足PH检验;模型整体检验P值0.416，模型整体满足PH检验。
test.ph <- cox.zph(res.cox);test.ph

##绘制Schoenfeld残差图以判断各协变量的PH假定是否成立。
#如图，实线是拟合的样条平滑曲线，虚线表示拟合曲线上下2个单位的标准差。如果曲线偏离2个单位的标准差则表示不满足比例风险假定。从图中可见，各协变量满足PH风险假设。
dev.new();ggcoxzph(test.ph)

##绘制Deviance残差图以判断Cox模型中的异常值。
#如图残差值均匀的分布在0上下，表明满足上述假定。
dev.new();ggcoxdiagnostics(res.cox,type = "deviance",
                 linear.predictions = FALSE,
                 ggtheme = theme_bw())
#观察各协变量残差分布
dev.new();ggcoxdiagnostics(res.cox,type = "dfbeta",
                 linear.predictions = FALSE,
                 ggtheme = theme_bw())


##绘制Martingale残差图以实现协变量的非线性诊断。
#如图，显示年龄局部有非线性趋势，但整体表现出线性趋势。
dev.new();ggcoxfunctional(Surv(time, status) ~ age+log(age)+sqrt(age),data = lung)




###————————————————————————示例1————————————————————————————
# Create the simplest test data set 
test1 <- list(time=c(4,3,1,1,2,2,3), 
              status=c(1,1,1,0,1,1,0), 
              x=c(0,2,1,1,1,0,0), 
              sex=c(0,0,0,0,1,1,1));head(test1)
# Fit a stratified model 
summary(coxph(Surv(time, status) ~ x + strata(sex), test1)) 

# Create a simple data set for a time-dependent model 
test2 <- data.frame(start=c(1,2,5,2,1,7,3,4,8,8), 
              stop=c(2,3,6,7,8,9,9,9,14,17), 
              event=c(1,1,1,1,1,1,1,0,0,0), 
              x=c(1,0,0,1,0,1,1,1,0,0));test2
summary(coxph(Surv(start, stop, event) ~ x, test2)) 

# Fit a stratified model, clustered on patients 
?bladder
bladder[1:25,]
bladder_1 <- bladder[bladder$enum < 5, ];head(bladder_1)
coxph(Surv(stop, event) ~ (rx + size + number) * strata(enum) + cluster(id), bladder_1)#enum即复发类型。本例中，一共有1234等不同的复发类型。探究复发类型与event的相关性时，医学上知道肿瘤复发与用药、肿瘤大小、肿瘤数目等有相关性，即enum与rx/size/number等协变量有多重共线性；另外，对于每种enum来说，rx/size/number的分布也不一定是均衡的。所以此处建模用的是(rx + size + number) * strata(enum)，即Cox回归的校正。

# Fit a time transform model using current age
lung[1:25,]
lung_cox <- coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung,tt=function(x,t,...) pspline(x + t/365.25))#这里自定义函数tt是可用于对时间变量进行变换。
summary(lung_cox)

tt <- function(x,t,...) pspline(x + t/365.25)
tt(75,365.25)
? pspline

###————————————————————————示例2(医学方)——————————————————————————————————————

library(coin)
data(glioma);head(glioma)

library(survival)
#the survival in subset of glioma 
g3 <- subset(glioma, histology =='Grade3')
#from different groups
fit <- survfit(Surv(time, event)~group,data = g3)
plot(fit, lty = c(2,1), col = c(2,1));legend('bottomright', legend = c('Control','Treatment'), lty = c(2,1), col = c(2,1))
#significance test;p= 0.0138
survdiff(Surv(time, event)~group,data = g3)
#Approximative Two-Sample Logrank Test:"|histology"指不同的组织类型中，分组的生存差别；以似然分布检验。
logrank_test(Surv(time, event)~group|histology,data = glioma, distribution = approximate(B = 1000))
str(glioma$histology)


data('GBSG2',package = 'TH.data')
head(GBSG2)#horTh指药物组，cens=1代表死亡
plot(survfit(Surv(time, cens)~horTh,data = GBSG2),lty = c(2,1), col = c(2,1), mark.time = T);legend('bottomright', legend = c('yes','no'), lty = c(2,1), col = c(2,1))
coxreg <- coxph(Surv(time,cens)~.,data = GBSG2);summary(coxreg)
#进一步挖掘数据。如图，发现：当pnodes
library(party)
tree <- ctree(Surv(time,cens)~.,data = GBSG2);dev.new();plot(tree)



###————————————————————————示例3:landmark分析————————————————————————————
library("landpred");library(landest);library(survival);library(survminer);library(ggplot2)

data(data_example_landpred);head(data_example_landpred,20)
data(example_obs);head(example_obs)
#TL:the observed event or censoring time for the primary outcome, equal to min(T, C) where T is the time of the primary outcome and C is the censoring time.
#DL:the indicator telling whether the individual was observed to have the event or was censored, equal to 1*(T<C) where T is the time of the primary outcome and C is the censoring time.
#TS:the observed event or censoring time for the intermediate event, equal to min(TS, C) where TS is the time of the intermediate event and C is the censoring time.
#DS:the indicator telling whether the individual was observed to have the intermediate event or was censored, equal to 1*(TS<C) where TS is the time of the primary outcome and C is the censoring time.
#Z:a baseline covariate vector
#treat:treatment indicator

## Our goal is to predict T𝕃 using information on T𝕊 and a discrete marker Z measured at baseline. For example, Z may be a binary marker indicating the presence of a gene mutation previously linked to the terminal event of interest.



###首先观察example_obs中不同的treat与生存结果的关系。从图中可以大致判断，当TL>2.5时，treat=1组有更低的累积生存率，而在此之前treat对累积生存率的影响几乎为0.
fit <- survfit(Surv(TL, DL) ~ treat, data = example_obs)
obs <- ggsurvplot(fit,
                           pval = TRUE, conf.int = F,
                           risk.table = F,
                           risk.table.col = "strata",
                           linetype = "strata",
                           surv.median.line = "hv",
                           palette = c("#E7B800", "#2E9FDF"))
p <- obs$plot + 
        theme(legend.justification=c(1,0),legend.position=c(0.95,0.85),
        axis.text.x = element_text(size = 10,colour = "black",face = "bold"),axis.text.y = element_text(size = 10,colour = "black",face = "bold"),#修改坐标轴标签属性
        axis.title.x = element_text(size = 15,colour = "black",face = "bold"),axis.title.y = element_text(size = 15,colour = "black",face = "bold")) + 
  scale_y_continuous(expand = c(0,0)) + 
  scale_x_continuous(expand = c(0,0)) #坐标轴原点重叠
dev.new();p
ggsave(p,file = 'example_obs.pdf')

##然而，在实际的过程中，是不是说treat是完全没有效果呢？此时要考虑intermediate event的发生。
delta.iptw.km(tl=example_obs$TL, dl = example_obs$DL, treat = example_obs$treat, tt=2,cov.for.ps = as.matrix(example_obs$Z))

delta.land.obs(tl=example_obs$TL, dl = example_obs$DL, treat = example_obs$treat, tt=2, landmark = 1, short = cbind(example_obs$TS,example_obs$DS), z.cov = as.matrix(example_obs$Z),cov.for.ps = as.matrix(example_obs$Z))#executable but takes time




#Estimates the probability of survival past some specified time using landmark estimation for an observational study setting
W.weight <- ps.wgt.fun(treat = example_obs$treat, cov.for.ps = as.matrix(example_obs$Z))
example_obs.treat <- example_obs[example_obs$treat == 1,]
surv.land.obs(tl=example_obs.treat$TL, dl = example_obs.treat$DL, tt=2, landmark = 1, short = cbind(example_obs.treat$TS,example_obs.treat$DS), z.cov = example_obs.treat$Z, ps.weights = W.weight[example_obs$treat == 1])

#设
置的这个是？
t0 = 2;tau = 8
#note: computationally intensive command below Prob.Covariate.ShortEvent(t0=t0,tau=tau,data=data_example_landpred)
out <- Prob.Covariate.ShortEvent(t0=t0,tau=tau,data=data_example_landpred)
out$data
data.plot <- out$data
dev.new();plot(data.plot$XS[data.plot$Z ==1], data.plot$Probability[data.plot$Z ==1], pch = 20, xlim = c(0,t0));points(data.plot$XS[data.plot$Z ==0], data.plot$Probability[data.plot$Z ==0],pch = 20, col = 2)


newdata <- matrix(c(1,1,0.5,1,0,3,0,1,1,1,4,1,1.5,1,0,10,1,5,1,0,11,0,11,0,1), ncol = 5, byrow=TRUE);newdata
#note: computationally intensive command below
out <- Prob.Covariate.ShortEvent(t0=t0,tau=tau,data=data_example_landpred,newdata=newdata)
out$newdata




###—————————————————————Cox回归相关函数————————————————————
#cluster()：如果每一组有多个成员，此时cluster可用于识别correlated groups of observations。如组别、窝别。
head(rats)
marginal.model <- coxph(Surv(time, status) ~ rx + cluster(litter), rats,subset=(sex=='f'));summary(marginal.model)
frailty.model <- coxph(Surv(time, status) ~ rx + frailty(litter), rats,subset=(sex=='f'));summary(frailty.model)

levels(cluster(rats$litter))
levels(frailty(rats$litter))


?rats#litter代表窝别

?frailty# add a simple random effects term to a Cox or survreg model
# Random institutional effect
coxph(Surv(time, status) ~ age + frailty(inst, df=4), lung)
coxph(Surv(time, status) ~ age + inst, lung)

# Litter effects for the rats data
rfit2a <- survreg(Surv(time, status) ~ rx + frailty.gaussian(litter, df=13, sparse=FALSE), rats,subset= (sex=='f'));rfit2a
rfit2b <- survreg(Surv(time, status) ~ rx + frailty.gaussian(litter, df=13, sparse=TRUE), rats,subset= (sex=='f'));rfit2b


#strata(): identifies stratification variables。比如不同的治疗这一分层变量。
?strata()
?ovarian
head(ovarian)
a <- factor(rep(1:3,4), labels=c("low", "medium", "high"));a
b <- factor(rep(1:4,3));b
levels(strata(b))
levels(strata(a,b,shortlabel=TRUE))
coxph(Surv(futime, fustat) ~ age + strata(rx), data=ovarian)#rx代表两种不同的治疗组。
ovarian$rx;strata(ovarian$rx)



###————————————————————Cox回归教学————————————————————————
library(survival)
library(splines)

#观察flchain的组成。flchain是研究游离轻链flc与生存时间的关系，有年龄、性别、取样年份、kappa、lambda等项。
set.seed(2018);flchain[sample(1:900,25),]


##一、不同性别对flc值的影响，在不同的年龄段是不一样的。即age会影响flc~sex。
#将年龄变换成分类变量：
age2 <- cut(flchain$age, c(49, 59, 69, 79, 89, 120),
              labels=c("50-59", "60-69", "70-79", "80-89", "90+"));age2
#计算flchain中每类age的数目：
counts <- with(flchain, table(sex, age2));counts
#自定义新列flc:
flchain$flc <- flchain$kappa + flchain$lambda
male <- (flchain$sex=='M')
mlow <- with(flchain[male,], smooth.spline(age, flc))
flow <- with(flchain[!male,], smooth.spline(age, flc))
#绘制年龄与flc含量的关系图：
dev.new();plot(flow, type='l', ylim=range(flow$y, mlow$y),
       xlab="Age", ylab="FLC");lines(mlow, col=2)
cellmean <- with(flchain, tapply(flc, list(sex, age2), mean, na.rm=T))#计算不同年龄、性别分组的flc的平均值。
matpoints(c(55,65,75, 85, 95), t(cellmean), pch='fm', col=1:2)#标注f/m
round(cellmean, 2)#保留2位小数输出。


###观察每个
options(contrasts= c("contr.treatment", "contr.poly")) # R default
cfit0 <- coxph(Surv(futime, death) ~ interaction(sex, age2), flchain)
cmean <- matrix(c(0, coef(cfit0)), nrow=2)
cmean <- rbind(cmean, cmean[2,] - cmean[1,])
dimnames(cmean) <- list(c("F", "M", "M/F ratio"), dimnames(counts)[[2]])#这里的list的函数是？
signif(exp(cmean),3)#计算e的cmean次方，结果保留3位数。


options(contrasts=c("contr.SAS", "contr.poly"))
cfit1 <- coxph(Surv(futime, death) ~ sex, flchain)
cfit2 <- coxph(Surv(futime, death) ~ age2 + sex, flchain)
cfit3 <- coxph(Surv(futime, death) ~ sex + strata(age2), flchain)#只要
# Unadjusted
summary(cfit1)
summary(cfit2)
anova(cfit2)
# Stratified
anova(cfit3)
summary(cfit3)
# Wald test
signif(summary(cfit2)$coefficients, 3)
anova(cfit1, cfit2)


###
casewt <- array(1, dim=c(2,5,4)) # case weights by sex, age group, estimator
csum <- colSums(counts)
casewt[,,2] <- counts[2:1,] / rep(csum, each=2)
casewt[,,3] <- rep(csum, each=2)/counts
casewt[,,4] <- 1/counts

wtindx <- with(flchain, tapply(death, list(sex, age2)))
cfitpop <- coxph(Surv(futime, death) ~ sex, flchain,
                   robust=TRUE, weight = (casewt[,,3])[wtindx])
cfityates <- coxph(Surv(futime, death) ~ sex, flchain,
                     robust=TRUE, weight = (casewt[,,4])[wtindx])
# Glue it into a table for viewing
tfun <- function(fit, indx=1) {
    c(fit$coef[indx], sqrt(fit$var[indx,indx]))
  }#？
coxp <- rbind(tfun(cfit1), tfun(cfit2,5), tfun(cfitpop), tfun(cfityates))
dimnames(coxp) <- list(c("Unadjusted", "Additive", "Empirical Population","Uniform Population"),c("Effect", "se(effect)"))
signif(coxp,3)


cfit4 <- coxph(Surv(futime, death) ~ sex * age2, flchain)
# Uniform population contrast
ysex <- c(0,-1, 0,0,0,0, -.2,-.2,-.2,-.2) #Yates for sex, SAS coding
contrast <- function(cmat, fit) {
  varmat <- vcov(fit)
  if (class(fit) == "lm") sigma2 <- summary(fit)$sigma^2
  else sigma2 <- 1 # for the Cox model case
  beta <- coef(fit)
  if (!is.matrix(cmat)) cmat <- matrix(cmat, nrow=1)
  if (ncol(cmat) != length(beta)) stop("wrong dimension for contrast")
  estimate <- drop(cmat %*% beta) #vector of contrasts
  ss <- qform(estimate, cmat %*% varmat %*% t(cmat)) *sigma2
  list(estimate=estimate, ss=ss, var=drop(cmat %*% varmat %*% t(cmat)))
}
contrast(ysex[-1], cfit4)























